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Abstract. The role of matter circulation between the disk and halo in establishing the volume filling factors of 
the different ISM phases in the Galactic disk {\z\ < 250 pc) is investigated, using a modified version of the three- 
dimensional supernova-driven ISM model of Avillez (2000). We carried out adaptive mesh refinement simulations 
of the ISM with five supernova rates (in units of the Galactic value), a/acai = 1, 2, 4, 8 and 16 (corresponding to 
starburst conditions) using three finer level resolutions of 2.5, 1.25 and 0.625 pc, allowing us to understand how 
resolution would affect the volumes of gas phases in pressure equilibrium. We find that the volume filling factors 
of the different ISM phases depend sensitively on the existence of a duty cycle between the disk and halo acting 
as a pressure release valve for the hot (T > 10^'^ K) phase in the disk. The mean occupation fraction of the hot 
phase varies from about 17% for the Galactic SN rate to ~ 28%, for a/aoai = 4, and to 44% for a/aoai = 16. 
The amount of cold gas (defined as the gas with T < 10"^ K) picked up in the simulations varies from a value of 
19% for a/aoai = 1 to ~ 5% for a/aoai = 4 and < 1% for higher SN rates. Background heating prevents the cold 
gas from immediate collapse and thus ensures the stability of the cold gas phase. The warm phase has volume 
fiUing factors varying between 0.25 and 0.37 for the three lowest SN rates used in the simulations. Overall the 
filling factor of the hot gas does not increase much as we move towards higher SN rates, following a power law of 
(/^ hot) ^ (""/""Gai)"'^^^. Such a modest dependence on the SN rate is a consequence of the evacuation of the 
hot phase into the halo through the duty cycle. This leads to volume filling factors of the hot phase considerably 
smaller than those predicted in the three-phase model of McKee & Ostriker (1977) even in the absence of magnetic 
fields. 
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1. Introduction 

In their seminal paper of a three-phase model regulated 
by supernova explosions in an inhomogeneous medium, 
McKee & Ostriker (1977) predicted a volume filling factor 
of the hot intercloud medium (HIM) of fy^hot — 0.7 — 0.8. 
However, observations point to a value of ^ 0.5 (e.g., 
Dettmar 1992) or even lower when external galaxies are 
taken into account (e.g.. Brinks & Bajaja 1986). A way 
out has been suggested by Norman & Ikeuchi (1989) by 
the so-called chimney model, in which hot gas can es- 
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cape into the halo through tunnels excavated by clustered 
supernova (SN) explosions. Indeed X-ray observations of 
several nearby edge-on galaxies have revealed extended, 
galaxy-sized halos (e.g., Wang et al. 2001). The transport 
of gas into the halo is, however, still controversial, and 
arguments that superbubble (SB) break-out may be in- 
hibited by a large-scale disk parallel magnetic field have 
been put forward by several authors (e.g., Mineshige et al. 
1993). 

It is true that magnetic tension forces can in prin- 
ciple considerably decelerate an upward flow, so that it 
will eventually stall before break-out. Tomisaka (1998) 
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has performed 3D MHD simulations of expanding super- 
bubbles including radiative cooling. He finds that bubble 
confinement only occurs when the energy injection rate is 
below a critical value of ~ lO'^^ergs"^ (see also MacLow 
& McCray 1988) and/or the scale height of the field is 
infinite. Coupling the field distribution to the vertical gas 
density distribution, as seems more realistic for frozen-in 
field lines, e.g. B oc p^^^, break-out will occur even if the 
midplane field strength is as high as 

Attempts to determine the occupation fraction of the 
different phases, and in particular of the hot gas, by means 
of modelling the effects of SNe and SBs in the ISM, have 
been carried out by several authors (e.g., Fcrricre 1995, 
1998; Korpi et al. 1999). However, these models do not 
include the circulation of gas between the disk and the 
full halo, thus being unable to resolve the high-z region; 
neither do they take into account the mixing between the 
different phases. Therefore, an estimate of the volume fill- 
ing factors may be misleading. 

Large scale simulations of the ISM in 2D (Rosen & 
Bregman 1995) and in 3D (Avillez 2000; Avillez & Berry 
2001) show the importance of the disk-halo interaction in 
the dynamics and evolution of the ISM in disk galaxies, 
which is intimately related to the vertical structure of the 
thick gas disk and to the rate of occurrence of supernovae 
per unit area in the Galactic disk. A major consequence 
of this interaction is the establishment of a duty cycle be- 
tween the disk and halo by hot gas breaking through the 
thick gas disk, cither in a violent way through chimneys 
or in a secular fashion through the buoyant rising of hot 
gas (Avillez 2000; Avillez & Mac Low 2001), into the up- 
per parts of the thick gas disk (the disk- halo interface). 
From here it escapes into the halo setting up a Galactic 
fountain (some of this gas may escape from the Galaxy as 
a wind, cf. Breitschwerdt et al. 1993) whose major frac- 
tion returns to the disk as cold gas (after condensing into 
clouds), interacting with the thin gas disk. 

In this paper we investigate how the volume filling fac- 
tors of the different ISM phases in the Galactic disk (which 
we consider hereafter to be the region with —250 < z < 
250 pc) evolve with the local Galactic SN rate (varying 
from the present value to a factor of 16 times higher) over 
a sufficiently long time scale in order to set up a dynami- 
cal cqiiilibrium condition. We show that matter circulation 
including the Galactic halo is an important ingredient in 
explaining the volume filling factors of the stable phases 
in the ISM, in particular a much lower value for the hot 
intercloud medium (HIM), even without magnetic fields, 
is obtained than advocated by McKee fc Ostriker (1977). 
The present study makes use of a modified version of the 
supernova-driven ISM model of Avillez (2000) and carries 
out simulations on grids of kpc-scale size with a spatial 
resolution as fine as 1.25 pc and 0.625 pc, respectively, in 
order to study the ISM in the disk and halo of star forming 
galaxies within a wide range of different supernova rates. 
The simulation time scales arc of the order of 400 Myr so 
that the disk-halo-disk cycle can be fully established, thus 



allowing the system to completely lose its memory of the 
initial conditions. 

In §2 we review the three-dimensional SN-driven ISM 
model used in the present work and the numerical schemes 
applied to it, as well as the set-up of the current study. 
In §3, we describe the global evolution of the ISM as seen 
in these runs. Section 4 discusses the results and their 
implications. Finally §5 presents a summary of the main 
results and final remarks. 

2. SN-Driven ISM modelling 

2.1. Numerics 

In the current work we report on large-scale simulations of 

the ISM including the disk and Galactic halo using a mod- 
ified version of the SN driven ISM model of Avillez (2000) 
coupled to a three-dimensional HD code using adap- 
tive mesh refinement (AMR) in a block-based structure, 
that relies on virtual topologies of CPUs created through 
Message Passing Interface (MPI) calls. This topology is 
associated to the computational domain that is divided 
into blocks, each having N=nxXn.yXnz cells. When a re- 
finement is required, a block is split into eight new blocks 
(children) of N cells. This process is repeated until the 
finest level of resolution is reached. All the information 
relative to the tree structure is preserved in the virtual 
topology, being only necessary to query the different CPUs 
to learn their location in this topology, and therefore, the 
location of their neighbors, children and parents. Each grid 
block is refined periodically in regions where steep pres- 
sure gradients appear. The local increase of the number 
of cells corresponds to an increase in linear resolution by 
a factor of two. The adaptive mesh refinement scheme is 
based on Bcrgcr & Colella (1989), but the grid genera- 
tion procedure follows that described in Bell et al. (1994). 
The gas dynamics part of the code uses the piecewise- 
parabolic method of Colella & Woodward (1984), a third- 
order scheme based on a Godunov method implemented 
in a dimensionally-split manner (Strang 1968) that relies 
on solutions of the Riemann problem in each zone rather 
than on artificial viscosity to follow shocks. 

2.2. Model 

2.2.1. Setup of basic processes 

The present model introduces substantial improvements 
to that of Avillez (2000), namely: (i) inclusion of back- 
ground heating due to starlight varying with z (Wolfire et 
al. 1995) and kept constant in the directions parallel to 
the Galactic plane {z = 0) where it is chosen to initially 
balance radiative cooling at 9000 K. With the inclusion 
of background heating the gas at T < 10^ K becomes 
thermally bistable; (ii) use of a tabulated cooling function 
taken from Dalgarno & McCray (1972) with an ionization 
fraction of 0.1 at temperatures below 10** K and a temper- 
ature cutoff at 10 K; (iii) inclusion of SNe type la with a 
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scale height of 325 pc (Freeman 1987); (iv) OB stars form 
in regions with a density and temperature thresholds of 10 
cm~^ and 100 K, respectively, and their mimbcr, masses 
and main sequence life times tms are determined from 
the initial mass function. Roughly sixty percent of the 
OB stars explode within the cluster, while the remain- 
ing (composed of the lowest mass stars whose location 
is determined kinematically by attributing to each star a 
random direction and a velocity) explode in the field with 
a scale height of 90 pc. The associations form in a layer 
with a scale height of 46 pc. The time interval between 
the explosions of all OB stars is determined by their tms 
for the given SN rate. 

As in the Avillez (2000), the present model includes 
the fixed vertical gravitational field provided by the stars 
in the disk and the interstellar gas is initially setup with 
a density stratification distribution that includes the cold, 
cool, warm, ionized and hot gases as described in Ferriere 
(1998). This ignores, by the way, curvature terms in the 
perpendicular r — z-plane, that may become relevant at 
some distance from the plane and induce radial motions. 
On the other hand, this effect will be small as long as 
the flow remains supersonic, since the gas has virtually no 
time to establish lateral pressure equilibrium. The model 
allows for the motion of the associations and includes mag- 
netic fields and cosmic rays. The latter two are the subject 
of forthcoming papers. 

2.2.2. Further processes 

Owing to the observational complexity of the ISM and the 

nonlinearity of the dominant processes occurring there, we 
have set up a model with the most important ingredients, 
neglecting other processes that may be relevant too, like 
differential rotation in the case of extended horizontal gas 
flows, self-gravity in case of Jeans unstable clouds, the 
gravitational field of a dark matter halo component, ther- 
mal conduction, and, most importantly, magnetic fields. 
The latter is the subject of a forthcoming paper (Avillez & 
Breitschwerdt 2004). Self-gravity is neglected, since with 
the cooling function we are using, we are not able to treat 
the formation of molecular clouds (also dust cooling is not 
included). In this respect our calculations are incomplete. 
However, it is known that the molecular phase decouples 
due to gravitational instability, and hence is not even in 
rough pressure equilibrium with the other ISM phases. 

The reason for attacking the problem in this way, in 
contrast to other authors who try to implement as many 
processes as possible with computing power being the ma- 
jor limitation, is the following. We believe that in order 
to gain insight into the nonlinear physics of a complex 
system like the ISM, driven by temporally and spatially 
variable energy input, we should first try to understand 
the interplay and effects of the most dominant processes. 
Therefore we do not strive for completeness here, but for a 
better understanding of the ISM gas dynamics. As a brief 
comparison with the observations shows, and, as will be 



shown in more detail in the second paper, the quantitative 
values that we obtain for the volume filling factors of the 
ISM gas are already in remarkably good agreement with 
the data. Below we discuss in some detail the relevance of 
processes that are known to be important in the ISM, but 
have been neglected in the present simulations. 

— Thermal conduction: The importance of thermal con- 
duction in the ISM has been the subject of debate for 
at least 25 years. There is some agreement that Spitzer 
conductivity does certainly overestimate the heat flux 
by a large amount, since the inclusion of an even wc;ak 
magnetic field reduces the mean free path of the elec- 
trons to a gyroradius, which is orders of magnitude 
below any astrophysical length scale. Saturated con- 
duction (e.g. Cowie & McKee 1977) is an improvement, 
when steep temperature gradients occur. In this case 
the heat flux will be independent of the temperature 
gradient. But the flux would be reduced by cos 9, where 
6 is the angle between the local magnetic field direc- 
tion and the temperature gradient. It has been argued 
recently in the context of cooling flows in clusters of 
galaxies (e.g. Narayan & Medvedev 2001), that turbu- 
lence may boost the conduction rate and thus over- 
come the suppression by a magnetic field. However, as 
has been pointed out by Malyshkin & Kulsrud (2001), 
there are two cifFects that inhibit heat flow even in the 
case of heat conduction in stochastic magnetic fields: 
(1) the electrons have to spiral along the lines of force, 
which are highly tangled, and therefore have to drift a 
long way, encountering on average smaller temperature 
gradients, (2) the electrons may become temporarily 
trapped on their journey by magnetic mirrors, until 
released by decreasing their pitch angle sufficiently by 
collisions. 

We have not included heat conduction in the present 
paper, since as a second order process it should be 
slower than hydro dynamic mixing. As is shown in this 
paper, supernova induced turbulence is quite strong in 
the ISM, and therfore we think that turbulent diffu- 
sion, induced by turbulent motions, will be the more 
efficient mixing process. 

— Coriolis forces: The present set-up ignores Galactic ro- 
tation, which can give rise to radial motions. A mea- 
sure of the importance of Coriolis forces with respect to 
inertial forces is the Rossby number cr = v/{Q.qL) 4C 
1, where L ~ 1 kpc is the maximum size of our grid. 
As far as the hot medium is concerned the rms veloc- 
ity component is of the order of 100 km/s and hence 
cr ~ 4, for fio = 8.4 x 10"-^^ s~-^. Moreover, the hot gas 
breaking out of the disk following the density gradient 
has a dominant z-component, for which the Coriolis 
term vanishes. 

— Dark matter halo: Up to a distance of ^ 10 kpc 
the amount of dark matter pulling the gas down is 
small compared to the disk potential. To see this, 
one can look at the gravitational potential, for which 
in case of the Milky Way one can use for the bulge 
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and disk component the mass model of Miyamoto & 
Nagai (1975) plus a spherical dark matter halo compo- 
nent as suggested by Innancn (1973). Comparing the 
gravitational acceleration terms near the solar circle, 
ffdisk and (7haio of the bulge plus disk and the dark 
matter halo, respectively, for a fluid element at ver- 
tical distances z = 5 and z = 10 kpc, one obtains 
[<7disk(^ = 5) = 1.19 X 10-8, ghaio(^ = 5) = 1.23 x IQ-^] 
and [5disk(^ - 10) = 8.72 x lO^^, (7haio(^ = 10) = 
1.96 x 10"^], in CGS units respectively. Only at z w 35 
kpc do the two terms become comparable, with the 
halo component dominating at larger heights. Thus ne- 
glecting the dark matter component in the simulations, 
where z < 10 kpc, is fully justified for spiral galaxies. 

2.3. Simulations 

The simulations use a computational domain that con- 
tains a section of the Galaxy with an area of 1 kpc^ and a 
vertical extension from -10 to 10 kpc. The innermost edge 
lies 8.5 kpc from the Galactic centre, coinciding with the 
solar circle. The grid resolution is 10 pc except in the layer 
between -500 and 500 pc, where two, three and four levels 
of AMR are used, yielding a finest level resolution of 2.5, 
1.25 and 0.625 pc, respectively. For the 2.5 and 1.25 pc 
resolutions five runs using cr/aQai = 1, 2, 4, 8 and 16 were 
carried out, while the 0.625 pc runs used a/acai = 1, 
2 and 4. At the cap surfaces z = — 10 kpc and z = 10 
kpc outflow boundary conditions were imposed and peri- 
odic ones along the four vertical boundary surfaces. The 
simulations were evolved for 400 Myr. Such a large grid 
size and long evolution time scale is needed to follow the 
time-dependent evolution of the ISM allowing several gen- 
erations of young massive stars to occur within the com- 
putational domain. The highest resolution simulations are 
used to answer the crucial question: does each refinement 
of the grid lead to a qualitatively and quantitatively differ- 
ent picture of the ISM or does the distribution of gas over 
the various phases converge to values depending mainly 
on the SN rate and not on grid resolution? 

3. Results 

3.1. Global evolution of the ISM 

The simulations start with a set of hydrodynamic and 
thermodynamic variables p, Pgas and T taken from obser- 
vations. However, the initially stratified distribution does 
not hold for long as a result of the lack of equilibrium be- 
tween gravity and (thermal, kinetic and turbulent) pres- 
sure during the "switch-on" phase of SN activity. As a 
consequence the gas distributions in the upper (z = 10 
kpc) and bottom (z = —10 kpc) parts of the grid col- 
lapse into the midplane, leaving low density material be- 
hind. When enough supernovae have gone off in the disk 
building up the required pressure support most of the col- 
lapsed gas expands supersonically through the grid filling 
all the computational domain. A continuous flow between 



the disk and halo is then set up, where the upward and 
downward flowing gas come into some sort of dynamical 
equilibrium. 

The gas in the disk also reaches equilibrium, although 
on a different time scale, mainly determined by the in- 
put of energy into the ISM by SNe, diffuse heating and 
the energy lost by adiabatic (due to expansion of SNRs, 
SBs and escaping disk gas) and radiative cooling. The disk 
equilibrium is related to the time scale that pressure needs 
to build up there as a result of SN explosions, while on 
the global halo scale such an equilibrium is only possi- 
ble after the full establishment of the duty cycle of the 
warm and hot gas and its circulation between the disk 
and the halo, which takes several hundred Myr. An up- 
per limit for this can be estimated by calculating the 
flow time, t/, that the gas needs to travel to the criti- 
cal point of the flow in a steady-state (see Kahn 1981). 
This is the characteristic distance from which information 
in a thermally driven flow can be communicated back to 
the sources. Then Tf ~ Tc/Cs, where Vc and Cg are the 
location of the critical point and the speed of sound, re- 
spectively. For spherical geometry, the critical point can 
be simply obtained from the steady state fluid equations, 
Tc ^ G Algal /i^c^), which yields with a Milky Way mass 
of Mgai « 4 X 10^^ M0 a distance Tc w 31.3 kpc for an 
isothermal gas at T = 2 x 10^ K (corresponding to a 
sound speed of 1.67 x lO^cms^^), and thus Tf ^ 180 Myr 
as an upper limit for the flow time. For comparison, the 
radiative cooling time of the gas at a typical density of 
n = 2 x 10"^ cm"^ is roughly Tc ^ 3kBT/{nA) « 155 Myr 
for a standard coUisional ionization equilibrium cooling 
function A = 8.5 x 10^^^ erg cm^ s~^ of gas with cosmic 
abundances. This value is apparently of the same order as 
the flow time, ensuring that the flow will not only cool by 
adiabatic expansion, but also radiatively, thus giving rise 
to the fountain return flow, which is the part of the out- 
flow that loses pressure support from below and therefore 
cannot escape. Note that Tc is the maximum extension 
of the fountain, more than an order of magnitude larger, 
than in simulations with grids with a vertical extension of 
1 kpc above and below the midplane. Those calculations 
definitely miss an important component of the galactic gas 
dynamics. 

In the disk, on the other hand, the cooling time scale 
is much shorter than in the halo, because of the higher 
ISM densities there. We thus find from our simulations 
that typical time scales of pressure fluctuations about a 
mean value of P/ks ^ 2600cm^'^K are of the order of 
30 Myr, which is by the way also a typical time scale for 
superbubble evolution. 

However, it should be emphasized, that since disk and 
halo are coupled dynamically not only by the escape of 
hot gas, but also by the fountain return flow striking the 
disk, the disk equilibrium will suffer secular variations on 
time scales of the order of 100 — 200 Myr. 

Cuts through the 3D data cube at z = pc and taken 
at 400 Myr presented in Figure Q] show the density (left 
column), P/k (middle column) and temperature (right col- 
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Left panel: mavillezJigla, Middle panel: mavillez_figlb, Right panel: mavillezjiglc 



Left panel: mavillez_figld, Middle panel: mavillez_figle, Right panel: mavillezjiglf 



Left panel: mavillezjiglg, Middle panel: mavillezjiglh, Right panel: mavillezjigli 



Fig. 1. Two dimensional cuts, through the 3D data cube, showing the p/ po (first column), P/k (middle column) and 
T (right column) distribution in the Galactic plane for a/acal = 1 (first row), 2 (second row) and 4 (bottom row). 



umn) distributions in the Galactic plane at 400 Myr of 
evolution, well after the system has reached dynamical 
equilibrium. The supernova rates in units of the Galactic 
value are: 1, 2, and 4, and the resolution of the finest 
AMR level is 1.25 pc. These images show the structure 
and morphology of the different ISM phases defined here 
as: cold (T < 10^ K), cool (10^ < T < 10'' K), warm 
(lO'' <T < 10^-^ K) and hot (T > 10^-^ K). The mor- 



phology of the ISM and the volume filling factors of the 
various phases obviously vary with SN rate. Not surpris- 
ingly, the amount of cold gas reduces with increasing SN 
rate, while the amount of the warm and hot gas increases. 
The images also show the expansion of SNRs and the in- 
teraction of shock waves, which can be best observed in the 
pressure maps where the shocks arc represented by high 
pressure and thin structures in red. This is not seen in 
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the density or temperature maps. A comparison between 
these maps shows that most of these shocks are located at 
the interfaces between hot and intermediate temperature 
gas, i.e., 10^ < T < lO^-^ K. 

Inspection of the pressure distribution shows that the 
most abundant colour is yellow-green, implying an average 
P /k ^ 2600 Kcm~^ (see also Fig. El upper right panel) 
for a = (JGai, and about 10,000 for a = AaQai- This is 
in agreement with the observational fact that the cool, 
warm and hot phase are in rough pressure equilibrium. 
This should however not be confused with a local equi- 
librium. On the contrary, severe deviations from the equi- 
librium values can occur locally due to SN activity and 
thermal instabilities (see the red an blue regions in Fig.^ 
middle column). 

For the Galactic SN rate, once the system reaches its 
equilibrium configuration, a thin disk of cold gas forms in 
the Galactic plane, and, above and below, a thick inhomo- 
geneous gas disk forms (this vertical structure of the disk 
is similar to that found in Avillez 2000). The code does 
not explicitly follow ionization states, but we can trace 
cool gas with a temperature T < 10^ K and a scale height 
of 180 pc, and warm gas with lO'' < T < 10^ K and a scale 
height of 1 kpc. These distributions reproduce those de- 
scribed in Dickey & Lockman (1990) and Reynolds (1987), 
respectively. The thick gas disk is punctured by chimneys 
that result from superbubbles occurring on either side of 
the midplane at a height greater than 100 pc. As they 
grow, they elongate along the z-direction, owing to the 
local density stratification of the ISM. Chimneys in the 
simulation typically have widths of approximately 150- 
200 pc. They inject high temperature gas directly from 
the Galactic disk into the halo, breaking through the cool 
and warm layers that compose the thick disk. This hot gas 
then contributes to the galactic fountain. Thus, the upper 
parts of the thick warm disk form the disk-halo interface, 
where a large scale fountain is set up by hot ionized gas 
escaping in a turbulent convective fiow. 

For higher SN rates a similar z-structure of a thick gas 
disk is seen, but with a higher vertical extension. This is 
mainly due to the rate of SNe in the field rather than to the 
clustered SNe, which drive superbubbles breaking through 
the thick gas disk and injecting their matter directly into 
the halo. The evolution of the vertical structure of star- 
forming galaxies with different SN rates will be discussed 
in more detail in a forthcoming paper. 



mavillez_fig2a 
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3.2. Probability distribution functions 

A major consequence of the set-up of the disk-halo cycle 
by SN activity is that the system loses any recollection of 
its initial conditions as most of the disk gas has already 
travelled into the halo and came back to the disk. Thus, 
the effects of initial conditions are not present in the tem- 
perature probability distribution functions (PDFs) of the 
system for different SN rates. 



Fig. 2. Averaged volume-weighted temperature PDFs 
over the periods of 0-50 Myr (red) and 350-400 Myr 
(black) calculated using 51 snapshots taken at time inter- 
vals of 1 Myr. The supernova rates used in these models 
are; (r/acai — 1 (top panel), 4 (middle panel), and 8 (bot- 
tom panel). The resolution of the finest AMR level is 1.25 
pc. 
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mavillezJigSa mavillez_fig4a 



mavillez_fig3b mavillez_fig4b 



mavillez_fig3c mavillezjig4c 



Fig. 3. Averaged volume- weighted density PDFs over the Fig. 4. Averaged volume-weighted pressure PDFs over the 

period 350-400 Myr calculated using 51 snapshots taken at period 350-400 Myr calculated using 51 snapshots taken at 

time intervals of 1 Myr for the SN rates shown in Figure[21 time intervals of 1 Myr for the SN rates shown in FigureEl 

The finest level AMR resolution is 1.25 pc. The finest level AMR resolution is 1.25 pc. 
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Figure 121 compares volume weighted histograms of the 
temperature over the periods of 0-50 Myr (red) and 350- 
400 Myr (black) calculated using 50 snapshots taken at 
time intervals of 1 Myr for the runs with a/acai — 1, 4, 
and 8 and a finest AMR resolution of 1.25 pc. These PDFs 
indicate that for low ct, the temperature peak is at about 
2000 K, making the cold/warm HI gas the most abundant 
gas phase, consistent with the density PDFs (Figure |3J 
and also with observations. This peak is shifted towards 
T « 40000 K, for a = l&aoaU making the warm phase 
the most important one. At the same time the relative 
importance of the hot phase increases as well, and one is 
moving towards a bimodal distribution. It can be directly 
seen from Fig. |21 that a simulation time of only 50 Myrs 
is not sufficient to obtain this result, mainly because the 
effect of upwards transport is to establish a duty cycle, 
which takes of the order of a few hundred Myrs. In any 
of the cases shown in the figures the averaged PDFs for 
the initial 50 Myr have two pronounced peaks, one around 
8000 K and the other around 5 x 10*^ K. Note that with the 
increase of SN rate the PDFs of the first 50 Myrs suffer 
large variations, indicating that the loss of recollection of 
the initial conditions will occur earlier for the highest SN 
rates, suggesting that in a time-asymptotic sense the SN 
rate is the controlling parameter. The corresponding av- 
eraged volume-weighted PDFs of the density and pressure 
distribution in the Galactic disk over the period 350-400 
Myr are shown in FiguresOland^ It can be seen that there 
are three distinct peaks in the density PDFs correspond- 
ing to the three most abundant regimes (cool, warm and 
hot). The cold regime has a peak that decreases steeply 
with increase of SN rate. These peaks have similar pressure 
ranges. The pressure PDFs show that the range of the to- 
tal pressure (dashed lines) decreases slightly with increase 
of the SN rate. For a = aaai and <dN/N>= 10"^ ^ 
tal pressure spans three orders of magnitude from 10^^^ 
to 10~^^, while for higher SN rates of 8 and 16 times the 
Galactic value the total pressure spans three and two or- 
ders of magnitude, respectively. This is indicative of the 
large variation in the pressure distribution between the 
different temperature regimes, and suggests that there are 
no real phases, i.e. co-existing thermodynamic regimes 
with different density and temperature but in pressure 
equilibrium. 

3.3. Volume filling factors 

At (7 = UGaU the hot phase has a volume filling factor 
fv,hot ^ 0.17, comparable to that of the cold gas, and the 
warm/cool phase is dominating. Increasing to 4, 8 and 16 
times acai pushes the warm gas to take over in volume, 
while there is a slight increase of the occupation fraction 
of the hot gas. It is interesting to note that the peak of the 
warm gas still surpasses that of the hot gas in magnitude 
even for a factor of 16 (see Fig. EI). This must be due to the 
fact that the reservoir of cold gas is used up by increased 
Lyman continuum photon absorption and shock heating. 



whereas the hot gas escapes into the halo. Between cr = 4 
and IQuGai, the cold/cool phase is reduced substantially 
and eventually to insignificance. This must have direct 
consequences for the formation of molecular clouds and 
thus for continuous star formation. It seems plausible that 
this leads eventually to self-regulation. 

Figure shows the history of of the volume filling 
factors of the different phases in the simulated disk for 
<^/<^Gai = 1,4, and 8, while Tableland Figure show the 
variation with supernova rate of the time averaged vol- 
ume filling factors of the different phases, over the period 
300-400 Myr calculated using 101 snapshots with a time 
interval of 1 Myr. 

The distribution of the occupation fraction of the dif- 
ferent phases comes to an equilibrium at about 100 Myr 
for ajoQai > 1, while for the Galactic SN rate it occurs 
somewhat later, i.e. at about 200 Myr. This is a result of 
the set-up of the disk-halo cycle, which for the Galactic 
SN rate takes about 200 Myr to establish, while for higher 
SN rates it takes less time in evacuating the hot gas from 
the disk. This is presumably due to an over-pressure with 
respect to the ambient medium and a higher sound speed 
of the injected hot gas. But in any case, the timescale for 
the cycle to be completed, even with higher SN rates, is 
always larger than 100 Myr, thus ruling out simulations 
with lower evolution times. 

Most remarkably, for ajaoai = 1, the hot (T > 10^-^ 
K) gas has a moderately low volume filling factor (it 
fluctuates around 0.17) in agreement with observations 
(^ 20%) even in the absence of magnetic fields, and is 
mainly distributed in an interconnected tunnel network, 
and in some cases it is even confined to isolated bubbles as 
seen in Figure ^ With the increase of the supernova rate 
to four and eight times the galactic rate the occupation 
volume of the hot gas increases to about 28% and 35%, 
respectively, after 400 Myr of disk evolution, which is still 
below the predictions of McKee & Ostriker (1977). Even 
for a/acai = 16, corresponding already to starburst con- 
ditions, the volume fraction of the hot gas increases only 
to 44% after 400 Myr of evolution. Such a behaviour can 
be approximated analytically by a power law fit (Figure [HI 

/ \ 0.363 

ifvMt) = 0.16 ^ , (1) 

with a a RMS percent error of 0.01. 

The warm gas becomes the dominant phase for 
(J /(jQai > 2. Between a joQai — 8 and 16 the occupation 
volume of this phase has a very small increase from 0.52 
to 0.54, as a result of the conversion of part of the cool gas 
into this phase. The decrease in the cool gas is a result of 
its conversion into the warm and hot phases. The growth 
of the hot phase in time is a mild one even with the large 
increase in the supernova rate. This is a consequence of 
the circulation of matter from the disk into the halo which 
acts as a pressure release valve for the disk gas. As the hot 
gas is vented into the halo, there is enough space for the 
10^ < r < 10^-^ K gas to be redistributed in the disk. 
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Table 1. Average volume filling factors of the different 
ISM phases for variable SN rate. The average was cal- 
culated using 101 snapshots (of the 1.25 resolution runs) 
between 300 and 400 Myr of system evolution with a time 
interval of 1 Myr. 



a" 


{fv,cold} 


{fvjCool} 




{fvjhot) 


1 


0.19 


0.39 


0.25 


0.17 


2 


0.16 


0.34 


0.31 


0.19 


4 


0.05 


0.30 


0.37 


0.28 


8 


0.01 


0.12 


0.52 


0.35 


16 


0.0 


0.02 


0.54 


0.44 



" SN rate in units of the Galactic SN rate. 
T < 10^ K. 
10^ < T < 10"* K. 
10" < T < 10^-^ K. 



T > 10^-^ K. 



mavillezJigSb 



mavillezJig6 



mavillezJig5c 



Fig. 6. Variation of the average volume filling factor of 

T < 103 (circles), 10^ < T < W (squares), 10^ < T < 
10^-^ (diamonds) and T > 10^-^ K (stars) phases in the 
simulated disk (for \z\ < 250 pc). The filling factors were 
averaged using 101 snapshots (of the 1.25 resolution runs) 
separated by 1 Myr between 300 and 400 Myr of the disk 
evolution. Red represents a power law fit to {fv,hot) with 
the power of 0.363. 



Fig. 5. Time cvohition of the volume filling factors of the 
cold (T < 10^ K), cool (10-'' < T < 10'' K), warm (10"* < 
T < 10^ -^ K), and hot (T > lO'"^-"^ K) phases for a/acai = 1 
(top panel), 4 (middle panel), and 8 (bottom panel). The 
finest AMR level resolution is 1.25 pc. Note that the line 
for the cold gas is slightly over = 0, for a/acai = 8. 



Furthermore, as there is recycling, low temperature gas 

continues to flow into the disk, thus contributing to the 
maintenance of the intermediate phases. 

For (J /(J Gal = 2 there seems to be a small difference 
in the volume occupation of the cool (aroimd 34%) and 
warm (around 31%) phases, while the occupation fraction 
of the cold gas differs from these by at most 15% of the 
total volume. It is interesting to note that the volume 
occupation of the cool gas decreases some 5% with the 
doubling of the SN rate for ujuGai < 4. On the other 
hand, there is a reduction of almost 70% of the filling 
factor of the cold gas when the supernova rate is increased 
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Fig. 7. Comparison between volume filling factors of cold 
(T < 10^ K), cool (10^ < T < 10"* K), warm (10'' <T < 
lO^-^ K) and hot (T > lO^-^ K) gas for the finest AMR 
level grid resolutions of 0.625 pc (black), 1.25 pc (red), and 
2.5 pc (green) for the SN rates a/acai = 2 (top panels), 
and 4 (bottom panels). 



to a/acai = 4. This fraction reduces to less than 1% at 
(j/<TQai = 8, disappearing completely for higher SN rates. 

3.4. Resolution effects 

Our simulations show how crucial spatial resolution is in 
order to capture small scale structures and promote the 
mixing of different fluid elements. For instance the amount 
of gas in the different temperature regimes, and most im- 
portantly in the cold phase, depends sensitively on spatial 
grid resolution. As an example consider the amount of 
cold gas, the volume filling factor of which increases for 
f/cGa; = 2 from ~ 9% in the 2.5 pc resolution calcu- 
lations to ^ 17% in the 1.25 pc and 0.625 pc resolution 



cases (s. Figure [7|l; and the discrepancy for the cold re- 
mains still large, when one increases the supernova rate to 
o'/'^Gai = 4. For the hot phase, on the other hand, there 
is not much difference in fv,hot with increasing resolution. 
The largest variation occurs for the smallest SN rates with 
a reduction of at most 20% compared to the values derived 
from the Ax = 2.5 pc resolution simulations. The resolu- 
tion increase also affects the cool and warm phases in a 
similar way with an increase in fv,cooi and a decrease in 

fv,warm fd' <^/<^Gal > 2. 

The variations in the volume filling factors are easily 
understood if one takes into account that with an increase 
of resolution it becomes possible to resolve the smallest 
scale structures. Thus instead of averaging out the gas 
density over larger cells, and thereby wiping out density 
peaks, radiative cooling as a nonlinear process can become 
more efficient, since high density regions contribute more 
to the energy loss rate than low density regions can com- 
pensate by an accordingly lower rate. Since cooling is most 
efficient for dense gas, the cool phase is affected most. In 
addition, the spatial resolution of shear layers and contact 
surfaces, gives rise to an increased level of turbulence and 
a larger number of mixing layers. The latter is most impor- 
tant, because it allows for a faster mixing between parcels 
of gas with different temperatures (conduction or diffusion 
processes being of second order and hence inherently slow 
in nature). The small scale mixing in these simulations is 
promoted by numerical rather than molecular diffusion, 
and therefore, the time scales for mixing in the different 
phases to occur is somewhat smaller (because it happens 
on larger scales) than those predicted by molecular diffu- 
sion theory (e.g., Avillez & Mac Low 2002). However, as 
already mentioned, turbulent diffusion, as a consequence 
of the onset of turbulence due to shear flows, will be most 
efficient. 

A comparison between the maximum density, Umax, 
and minimum temperature, Tmm, measured at the dif- 
ferent finer level resolutions reveals that an increase in 
resolution from Ax = 2.5 to 1.25 pc implies an average 
increase in Umax and a decrease in Tmin by factors greater 
than 5 at any SN rate (cf. FigurelHl. When resolution is in- 
creased from 1.25 to 0.625 pc the differences between Umax 
and Tmin for all the SN rates are small and diminish with 
the increase of SN rate. The high Tmin for cr/acai — 4 
at any resolution is due to the higher cooling time of the 
gas supply with on average higher temperature. Figure El 
compares the average values of Tmin and Umax calculated 
between 200 and 400 Myr of evolution for the three reso- 
lutions (for details, see figure caption). 

When a resolution of 0.625 pc is used the average mini- 
mum temperature and maximum density, Tmin and Umax , 
suffer small variations with respect to those determined 
at Ax = 1.25 pc simulations for all the SN rates dis- 
cussed above. This is a clear indication of the convergence 
of the simulations in reproducing the physical processes 
involved in the dynamics and evolution of the ISM and 
that the simulations with Ax < 1.25 pc can represent the 
real ISM. From the fit for (T™„) = 8.35 exp( j-^f^), and 
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Left panel: mavillezJigSa, Middle panel: mavillez_fig8b, Right panel: mavillezjig8c 



Fig. 8. Comparison between grid resolutions of 0.625 pc (black), 1.25 pc (red) and 2.5 pc (green) of maximum density 
(top panel) and minimum temperature (bottom panel) for a/acal = 1 (left column), 2 (middle column), and 4 (right 
column). Note that in any of the panels there is large difference between the 2.5 and 1.25 pc resolutions, while the 
differences between 0.625 pc and 1.25 pc resolutions are much smaller. This is an indication of the convergence of the 
simulations in reproducing the physical processes involved in the dynamics and evolution of the ISM, thus representing 
a fairly realistic distribution of gas in density-temperature space. 



(iimaa:) = 442 exp(- ^^f ^ ), respectively, we can see that 
there is rapid convergence for Aa; ^ 1.1 pc. 

4. Discussion 

The filling factors of some of the phases obtained in this 
study are similar to those estimated by Spitzer (1990), 
Ferriere (1995) and Avillez (2000), but seem to be in con- 
tradiction with predictions of Ferriere (1998) and are not 
consistent with Korpi et al. (1999). 

The first three authors estimated a volume filling fac- 
tor of the hot gas in the Galactic disk of ^ 0.2, while for 
the cold and neutral phases Ferriere (1995) and Avillez 
(2000) arrived to similar values. The largest variation be- 
tween the volume filling factors in Avillez (2000) and the 
present paper is by far in the value of fv,coid, as a result 
of the introduction of background heating due to starlight 
leading to the formation of thermally stable branches at 
103.9 <T< 104-2 K and T< 100 K. Most of the cold gas is 
located in the unstable branch at temperatures below 10^ 
K (Avillez & Breitschwerdt 2004). 

Based on models which follow the expansion and con- 
traction phases of SNRs and SBs, Ferriere (1998) esti- 
mated the variation of the volume filling factor of the hot 
gas with z. It is argued that near the solar circle, the 



volume occupation of the hot gas increases from 0.15 at 
2; = pc to 0.23 at z = 200 pc and decreases gradually for 
z > 200 pc (Figure 12 of Ferriere 1998). The discrepancy 
between these results and the calculations reported in this 
paper can be understood by the fact that Ferriere's cal- 
culations reflect the distribution of the hot gas inside the 
SNR and SB cavities whose number decreases with z, and 
therefore the volume filling factor of hot gas decreases ac- 
cordingly. This contrasts with the present paper where in 
a global 3D simulation the hot gas is not exclusively con- 
nected to individual SNRs or SBs^, but instead can rise 
into the halo. As the hot gas breaks out of the thick disk, 
in which its volume filling factor first decreases, it even- 
tually starts filling the entire volume of the Galactic halo. 
Thus, Ferriere's model does in fact not resolve the vertical 
gas transport, and hence misses out on the volume frac- 
tion of a considerable amount of gas, leading effectively to 
an increase with z-height above 1 kpc as described here. 

Figure 2 of Korpi et al. (1999) shows the temperature, 
density and P/k PDFs for different volumes of the disk: 
\z\ < 0.25, 1^1 < 0.5, and \z\ < 1 kpc. These are averaged 

^ In fact we see a large number of SNRs and SBs disintegrat- 
ing on timescales of 1 to 30 Myrs, respectively, even inside the 
Galactic disk, whore density gradients are on average smaller 
than in the halo. 
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Fig. 9. Comparison between the average minimum temperature (left panel) and maximum density (right panel) as 
function of resolution for the three SN rates: a/aGal = 1 (open circles), ujuGal = 2 (stars), and a/acal = 4 (triangles). 
The plots also show exponential fits to the data points. These average values were calculated during the last 200 Myr 
of evolution, i.e. after establishing dynamical equilibrium, so that their history does not reflect any memory of the 
initial conditions. 



volume-weighted histograms calculated over a period of 50 
Myr using six snapshots at equal time intervals between 
20 and 70 Myr of evolution. The temperature PDF for 
the Galactic disk gas has a bimodal distribution with two 
strong peaks one at 10^ K and another at ~ 10^ K. The 
latter corresponds to a fv,hot varying from 20 — 30% at 
2; = pc to 50 — 60% for z = 300 pc. These values are 
in disagreement with those discussed in the present paper 
for the Galactic SN rate. The discrepancy is related to 
the small grid extent in the ^-direction (with — 1 < ^ < 1 
kpc) used by these authors, and consequently the increas- 
ing escape of material from the grid with time without 
any return flow making their model, as the authors argue 
themselves, only meaningful for a limited length of time, 
which as we discuss below should not exceed 100 Myr. 



The need for a duty cycle and the establishment 
of a global dynamical equilibrium require the use of 
an extended grid in the direction perpendicular to the 
Galactic plane. The lack of such an extended 2;-grid in- 
hibits the disk-halo-disk circulation of matter, which oth- 
erwise would return gas to the disk sometime later, with 
noticeably increasing effects for the dynamical evolution 
as time proceeds, which can be clearly seen in the present 
simulations. In grids, which have a small vertical exten- 



sion, e.g., only 1 kpc above and below the midplane^, com- 
pared to the maximum height to which the hot plasma 
rises due to the injection of energy and momentum from 
the sources, most of the gas escapes from the disk in less 
than 100 Myr (without ever returning to it). Thus the 
simulations can only be followed for a small evolution 
time, and therefore, the duty cycle cannot be set up and 
the system never reaches a dynamical equilibrium state. 
Moreover, the volume weighted histograms of the thermo- 
dynamic properties (e.g., density, pressure and tempera- 
ture) of the disk gas will retain a memory of the initial 
evolution. As explained in §3.1, after an initial collapse of 
the matter distribution towards the grid midplane, pres- 
sure is built up and there is a redistribution of matter 
on the grid, filling it and allowing a substantial fraction 
of the gas to escape through the top and bottom bound- 
aries. Therefore, averaged PDFs that including 20 - 50 
Myrs snapshots will show the presence of a large fraction 
of cold gas from the collapse phase (recall the averaged 
50 Myr PDF shown in the top panel of Figure 2 in §3.2), 
while PDFs constructed at later times will show the dom- 
inance of the hot gas. A combination of these PDFs would 

Note, that the maximum extension of the fountain is more 
than an order of magnitude larger, and therefore such restricted 
calculations definitely miss an important component of the 
galactic gas dynamics. 
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give a pronounced bimodal distribution (even if there is 
no background heating due to starhght) corresponding to 
large volume fiUing factors for the hot and coolest ISM 
phases. In other words such PDFs are a reflection of the 
initial evolution of the simulations and not of the time 
when a dynamical equilibrium is already established. 

5. Summary and final remarks 

In the present paper we have described a set of 3D hydro- 
dynamical simulations of the interstellar medium in order 
to study the distribution of the ISM phases and how they 
vary with increasing star formation rate (i.e., SN rate) in 
star forming galaxies. The major goals of this work were 
(i) to see if the presence of a disk-halo-disk circulation has 
a major impact on the volume filling factors of the hot 
phase, and (ii) which minimum grid resolution is needed in 
order to obtain quantitatively reliable results that can be 
compared to observations. The main results of the present 
work can be summarized as follows: 

— The occupation fractions of the different ISM phases 
depend sensitively on the presence of a duty cycle es- 
tablished between the disk and halo working as a pres- 
sure release valve for the hot phase. 

— The mean occupation fraction of the hot phase varies 
from about 17% for the Galactic SN rate to 28% and 
44% for (j/<JGai — 4, and 16, respectively. The mean 
occupation fraction follows a power law increase with 
SN rate, with the power law index of ^ 0.363. 

— The amount of cold gas picked up in the simulations 
varies from a value of roughly 19% (for the Galactic 
SN rate) to about 5% for a/aoai = 4 and < 1% for 
o-fa-Gai > 8. 

— The warm phase has occupation fractions varying be- 
tween 25% and 37% for the three lowest SN rates 
{<j/crGai < 4). 

— Background heating is the main reason for the increase 
in the amount of cold gas in comparison to that in 3D 
simulations without any background heating. 

— A minimum grid resolution of 1.25 pc is needed for 
quantitatively reliable results, as the convergence to- 
wards the 0.625 pc resolution simulations shows. 

— A SN rate of 8 — 16 times the Galactic value already 
represents a starburst; there is increasing evidence that 
most star forming galaxies have undergone several such 
phases during their evolution. In particular for high 
redshift galaxies, starbursts seem to have been com- 
mon. Our simulations show that the volume filling fac- 
tor for the hot phase in the disk increases only moder- 
ately from {fv,hot) — 0.35 and 0.44 for oja^^x = 8 and 
16, respectively. 

This implies that in X-ray observations, the value 
of {fv,hot) in the disk is not a reliable indicator for 
starburst. Instead the size of the halo in soft X-rays 
is strongly correlated with a starburst as can be seen 
from the size of the X-ray halos in recent XMM- 
Newton observations of NGC 253 (e.g., Pietsch et al. 



2001) and NGC 3079 (Breitschwerdt et al. 2004), and 
also ROSAT observations of M 82 (e.g., Bregman et 
al. 1995) and NGC 253 (Dahlem et al. 1998, Pietsch 
et al. 2000). 

— Even for Galactic SN rates the fountain cycle is es- 
tablished and thus hot gas is present in galactic ha- 
los. This explains why also star forming galaxies with 
SN rates comparable to the Galaxy exhibit soft X-ray 
halos, albeit smaller than in starburst galaxies, as has 
been observed for NGC 4631 (Wang et al. 2001), which 
also may have been disturbed by a companion galaxy 
thereby enhancing the star formation rate, and NGC 
891 (Bregman & Irwin 2002), which is often referred 
to as the twin galaxy to the Milky Way. 

The calculations presented in this paper do not include 
the magnetic field. A parameter study of the effects of the 
magnetic field and cosmic rays in the ISM is underway 
and will be described in Avillez & Breitschwerdt (2004). 
It will deal with the effects of the B-field dissipation, which 
is too low for Galactic dynamos (Ferriere 1998), as well as 
with the effects of a weak and strong magnetic field. If the 
magnetic field is present and is initially mainly oriented 
parallel to the disk, transport in the halo may be inhib- 
ited, although not prevented. On larger scales magnetic 
tension forces become weaker than on the smallest scales 
and therefore vertical expansion might still take place effi- 
ciently. Either way the occupation fraction of the hot gas 
could be comparable to the values observed in the present 
simulations. 
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